In this lab we will explore and implement some of the most common processes that are applied to images from raw to JPG or other formats. We will also use this to get familiar with some simple image processing in Python.
We will use scipy.ndimage to load and save images which have a simple interface and stores the data as simple numpy arrays.
You are free to change the code I give you or not to use it at all.
In [1]:
#some magicto show the images inside the notebook
%pylab inline
%matplotlib inline
import matplotlib.pyplot as plt
from scipy import ndimage
import numpy as np
# A hepler function for displaying images within the notebook.
# It displays an image, optionally applies zoom the image.
def show_image(img, zoom=1.5):
dpi = 77
plt.figure(figsize=(img.shape[0]*zoom/dpi,img.shape[0]*zoom/dpi))
if len(img.shape) == 2:
img = np.repeat(img[:,:,np.newaxis],3,2)
plt.imshow(img, interpolation='nearest')
# A hepler function for displaying images within the notebook.
# It may display multiple images side by side, optionally apply gamma transform, and zoom the image.
def show_images(imglist, zoom=1, needs_encoding=False):
if type(imglist) is not list:
imglist = [imglist]
n = len(imglist)
first_img = imglist[0]
dpi = 77 # pyplot default?
plt.figure(figsize=(first_img.shape[0]*zoom*n/dpi,first_img.shape[0]*zoom*n/dpi))
for i in range(0,n):
img = imglist[i]
plt.subplot(1,n,i + 1)
plt.tight_layout()
plt.axis('off')
if len(img.shape) == 2:
img = np.repeat(img[:,:,np.newaxis],3,2)
plt.imshow(img, interpolation='nearest')
lighthouse = "data/lighthouse_RAW_noisy_sigma0.01.png"
img=ndimage.imread(lighthouse)
show_image(img[200:400,200:400],3) #takes a region of the image
The first task is to get familiar with the tool DCRAW. www.cybercom.net/~dcoffin/dcraw/. You are given an image _MG_4257.CR2 which is the RAW format for Canon DSLR cameras. Install it in your computer and run it, you will get the options and usage. In this task you should define some functions that simply call the dcraw command with the appropriate options from the a python programm.
Convert to a normal jpg in sRGB color space. Test different options (for example GAMMA) and compare the results. You are given a JPG file for this image which is the one produced by the camera
Identify the option for this and recover the original raw data in a standard image format (pgm). This will be a black and white image with a bayer pattern.
The standard pipeline applies color and gamma corrections, when we use images for computation is best to have the linear input and preserve as much information as possible. Find the options to recover the linear 16 bits tiff image from the raw.
Capture an image with your camera in RAW mode, extract the RAW image using the command line software dcraw (don’t process the image, just extract the linear, RAW),
In the next tasks test your code for demosaicking and gamma correction on your own raw images, compare with the results you get from dcraw (for demosaicking and conversion to sRGB). If you can’t capture, you have an example in data.
In [2]:
###HINT: Look for the "subprocess" module
import subprocess, shlex
path="data/raw/_MG_4257.CR2"
def raw2PPM(path_to_file, g1,g2):
g1=str(g1)
g2=str(g2)
my_cmd = './dcraw -o 1 -g '+g1+' '+ g2 + ' -c '+path_to_file + ' > image1'
args = shlex.split(my_cmd)
subprocess.call(args)
subprocess.Popen(my_cmd, shell=True)
def extract_raw(path_to_file):
my_cmd = './dcraw -o 0 -D -c '+path_to_file + ' > image2'
args = shlex.split(my_cmd)
subprocess.call(args)
subprocess.Popen(my_cmd, shell=True)
def raw2linear_tiff(path_to_file):
my_cmd = './dcraw -4 -T -c '+path_to_file + ' > image3'
args = shlex.split(my_cmd)
subprocess.call(args)
subprocess.Popen(my_cmd, shell=True)
g1=4
g2=1
raw2PPM(path,g1,g2)
extract_raw(path)
raw2linear_tiff(path)
In [3]:
img1=ndimage.imread("image2")
#show_image(img[1000:1800,1300:2000],1) #takes a region of the image
show_image(img1[1000:1800 , 1300:2000],0.2) #takes a region of the image
Take the 16 bits linear image and apply a gamma encoding to it. https://en.wikipedia.org/wiki/Gamma_correction
Display the linear image, the gamma encoded one, and the JPG produce by dcraw.
In [4]:
# Verbose aliases for gamma transformation
def gamma_decode(img, gamma=2.2):
#Input is a gamma encoded image, output is a linear image
scaling=np.power(img.astype(float),np.ones(img.shape)*gamma)
normalizing=scaling / np.linalg.norm(scaling)
return normalizing*np.iinfo(img.dtype).max
def gamma_encode(img, gamma=2.2):
#Input is a linear image, output is a gamma encoded image
scaling=np.power(img.astype(float),np.ones(img.shape)/gamma)
normalizing=scaling / np.linalg.norm(scaling)
return normalizing*np.iinfo(img.dtype).max
test_encode=gamma_encode(tests[0])
test_decode=gamma_decode(tests[0])
show_images([test_encode,test_decode])
#load the converted images and display using show_images()
White balance is the process of setting whites to white. Depending on the lighting, white pixels will look of the colour of the light, in our example whites look quite yellowish. https://en.wikipedia.org/wiki/Color_balance
Different methods exist to estimate the color of the light, you will implement two simple white balance methods.
The white point simple takes in a pixel (three values) which we know is white, for example a pixel from the fence, and use this pixel to set it to white, and the rest of the image accordingly.
The second method is the "grey world assumtion". The idea here is to compensate the color channels so their average is grey (all channels have the same value). To do this, compute the mean of each channel (RGB), normalize to the maximum value of the triplet, and use this as your compensation factor.
In [5]:
def white_balance_white_point(img, whitePoint):
img=img/whitePoint
return img #!change accordingly
def white_balance_gray(img):
av=np.average(np.average(img, axis=0), axis=0)
maxa=np.amax(av)
alfa=av/maxa
img=img/alfa
return np.clip(img,0.,1.)
im = ndimage.imread("data/lighthouse.png")/255.
im1 = copy(im)
im2 = copy(im)
im_wp = white_balance_white_point(im1,im1[309:310,348:349,:].flatten())#im[309:310,348:349,:].flatten()) #you can select a different pixel
im_g = white_balance_gray(im2)
show_images([im,im_wp,im_g],0.5)
np.clip([ 0.29411765, 3.64705882 , 3.68627451],0.,1.)
Out[5]:
In [ ]:
#White balance other images and display using show_images()
Implement demosaicking using linear interpolation. The example image ‘lighthouse_RAW_colorcoded.png’ shows you in false colors which subpixel measures with color channel. Don’t use that for demosaicking, use the image ‘lighthouse_RAW_noisy_sigma0.01.png’ Check your implementation with other images from data/raw. Lecture 2 slides here http://franchomelendez.com/Uwr/teaching/COMPHO/_LECTURES/L2/digital_photography.html
Unfortunatelly, not all the sensors have the same pattern, so you will find that the patern for lighthouse image is not valid for the other raw images. In essence in these cases the result is that the blue and red channels are exchanged. You can ignore this for the naive solution and hardcode the pattern, but you should take it into account when developing the other solutions.
In [6]:
# A helper function for merging three single-channel images into an RGB image
def combine_channels(ch1, ch2, ch3):
return np.dstack((ch1, ch2, ch3))
img = ndimage.imread("data/lighthouse_RAW_noisy_sigma0.01.png")
def deX(ch3, start):
for i in range(start,len(ch3),2):
for j in range(start,len(ch3[0]),2):
sum=0
div=4
try:
sum+=ch3[i-1][j-1]
except IndexError:
div-=1
try:
sum+=ch3[i-1][j+1]
except IndexError:
div-=1
try:
sum+=ch3[i+1][j-1]
except IndexError:
div-=1
try:
sum+=ch3[i+1][j+1]
except IndexError:
div-=1
ch3[i][j]=sum/div
return ch3
def deH(ch2, start):
for i in range(0,len(ch2)):
for j in range((i+start)%2,len(ch2[0]),2):
sum=0
div=0;
if(i>1):
sum+=ch2[i-1][j]
div+=1
if(i+1<len(ch2)):
sum+=ch2[i+1][j]
div+=1
if(j>1):
sum+=ch2[i][j-1]
div+=1
if(j+1<len(ch2[0])):
sum+=ch2[i][j+1]
div+=1
ch2[i][j]=sum/div
return ch2
def demosaic_simple_naive(img):
ch1=img*1
ch2=img*1
ch3=img*1
ch2=deH(ch2,0)
#print(ch2-img)[:6,:6]
ch3=deX(ch3,0)
ch3=deH(ch3,1)
#print(ch3-img)[:6,:6]
ch1=deX(ch1,1)
ch1=deH(ch1,1)
return combine_channels(ch1,ch2,ch3) #!change accordingly
print img[:6,:6]
show_image(demosaic_simple_naive(img)[:,:]) #should show a color image
#print(combine_channels(img,img,img))
In [7]:
image=copy(img)
show_image(np.maximum(np.ones(np.shape(image)),(2**16-1)-image))
test_image=np.maximum(np.ones(np.shape(image)),(2**16-1)-image)
tests=[img,test_image]
So far you didn't need to think too much, this should be a bit more challenging. Some things that might be useful are numpy.tiles which replicates a pattern a given number of times.
You can use this for example to specify mask, which you can simply multiply by the image to get a masked image. Masks should have only values of 0 or 1.
Convolutions are an important concept in image processing. You have a visual explanation here http://setosa.io/ev/image-kernels/ convolutions are expensive operations, specially when the kernels get big. Fortunatelly ndimage provides us with implementations that work like this:
kernel = np.array([[1,0,1],[0,0,0],[1,0,1]])
#This kernel only uses the corner pixels
gradient = ndimage.convolve(image, kernel)/4
#In order to normalize the value of the pixel, we need to divide by the sum of the values in the kernel
#This would be equivalent to:
kernel = np.array([[0.25,0,0.25],[0,0,0],[0.25,0,0.25]])
gradient = ndimage.convolve(image, kernel)
You can use convolutions to compute the gradient in the next task as well.
In [10]:
patterns_bayer='RGGB'
def recognize_pattern(pattern,shape):
s=shape
i1,i2,i3,i4=(0,0,1,1)
if pattern=='BGGR':
i1,i2,i3,i4=(1,1,0,0)
elif pattern=='GBRG':
i1,i2,i3,i4=(1,0,0,1)
elif pattern=='GRBG':
i1,i2,i3,i4=(0,1,1,0)
v1=np.array([int(x%2==i1) for x in range(s[0])]).reshape(-1,1)
v1t=np.array([int(x%2==i2) for x in range(s[1])]).reshape(1,-1)
v2=np.array([int(x%2==i3) for x in range(s[0])]).reshape(-1,1)
v2t=np.array([int(x%2==i4) for x in range(s[1])]).reshape(1,-1)
#print(v1,v1t,v2,v2t)
red=v1*v1t
blue=v2*v2t
green=(np.ones(s)-red-blue).astype(uint8)
return red,green,blue
def demosaic_simple_rgb(img,pattern='RGGB'):
typ=img.dtype
s=img.shape
kek=np.ones(s)
v1m,gp,v2m=recognize_pattern(pattern,s)
#print(v1,v1t,v2,v2t)
#brdpatt=v1m+v2m
#gpatt=np.ones(s)-brdpatt
blue=img*v2m
red=v1m*img
green=gp*img
kernelrb=np.array([[0.25,0.5,0.25],[0.5,1.0,0.5],[0.25,0.5,0.25]])
kernelg=np.array([[0,0.25,0],[0.25,1.0,0.25],[0,0.25,0]])
red2=ndimage.convolve(red,kernelrb).astype(typ)
green2=ndimage.convolve(green,kernelg).astype(typ)
blue2=ndimage.convolve(blue,kernelrb).astype(typ)
return (red2,green2,blue2)
def demosaic_simple(img,pattern='RGGB'):
### Reimplement simple demosaicing without using loops
red,green,blue=demosaic_simple_rgb(img,pattern)
return combine_channels(red,green,blue) #!change accordingly
kek=demosaic_simple(tests[0])
print(kek.dtype)
show_image(demosaic_simple(tests[0])/255.0)
Implement Edge-based Demosaicing. Compare the results. Details are in the lecture notes.
For the green channel use this recipe:
Summary:
For the other two channels use green based interpolation (example for the red channel):
Do the same for blue
In [12]:
tests.append(tests[0][500:518,400:418])
def demosaic_green_edge_based(img, patterns=patterns_bayer, offset=(0,0)):
### Implement the green edge based demosaicing algorithm
typ=img.dtype
s=img.shape
rp,gp,bp=recognize_pattern(patterns,s)
greenbase=img*gp
greenh=abs(ndimage.convolve(greenbase.astype(int),[[1,0,-1]])).astype(typ)
greenv=abs(ndimage.convolve(greenbase.astype(int),[[1],[0],[-1]])).astype(typ)
gbool=greenh>greenv
gbool2=greenh<greenv
gbool3=greenh==greenv
green=ndimage.convolve(greenbase,[[0, 0, 0],[0.5, 1.0, 0.5],[0, 0, 0]])*gbool2#*(rp+bp)
green=green+ndimage.convolve(greenbase,[[0, 0.5, 0],[0, 1.0, 0],[0, 0.5, 0]])*(gbool)#*(rp+bp)
green=green+ndimage.convolve(greenbase,[[0, 0.25, 0],[0.25, 1.0, 0.25],[0, 0.25, 0]])*(gbool3)#*(rp+bp)
green=green.astype(typ)
red=((img).astype(int)-green)*rp
blue=((img).astype(int)-green)*bp
kernelrb=np.array([[0.25,0.5,0.25],[0.5,1.0,0.5],[0.25,0.5,0.25]])
red=np.clip((ndimage.convolve(red,kernelrb)+green),0,np.iinfo(img.dtype).max).astype(typ)
blue=np.clip((ndimage.convolve(blue,kernelrb)+green),0,np.iinfo(img.dtype).max).astype(typ)
return combine_channels(red,green,blue)
lh_basic = demosaic_simple(img)
lh_gb = gamma_encode(demosaic_green_edge_based(img))
#this should show the images demosaiced
show_images([lh_basic, lh_gb], zoom=0.5, needs_encoding=True)
show_images([lh_basic[450:550, 50:150],
lh_gb[450:550, 50:150]], zoom=3.5, needs_encoding=True)
Unfortunately, the lighthouse is a difficult example for demosaicking. You should see color artifacts in the fence. Go back to the linear demosaicked image, convert it to YCrCb color space and median filter the chrominance channels but not the luminance channel. Convert back to RGB and then apply the sRGB gamma. You can use the implementation of the median filter available in ndimage, but you need to understand how it works and how you could implement a simple version.
You can find the coefficients in JPEG specification (https://www.w3.org/Graphics/JPEG/jfif3.pdf).
The YCbCr color space convertion needs to add a constant (128 or 0.5) to avoid negative values, since image formats don't allow for negative values. However, in our case, we don't care if we have negative values in our represntation (since they are floats) and because we are converting back to RGB. So you can safely ignore those constants for this task, which simplifies the convertion to a dot product.
In [108]:
import scipy
# Functions for converting color space on a three-channel image between RGB and YCbCr.
#Y = 0.299 R + 0.587 G + 0.114 B
#Cb = - 0.1687 R - 0.3313 G + 0.5 B
#Cr = 0.5 R - 0.4187 G - 0.0813 B
#R = Y + 1.402 (Cr)
#G = Y - 0.34414 (Cb) - 0.71414 (Cr)
#B = Y + 1.772 (Cb)
def rgb2ycbcr3(dtuple):
R,G,B=dtuple
Y = 0.299* R + 0.587 *G + 0.114 *B
Cb = - 0.1687* R - 0.3313* G + 0.5 *B #+128
Cr = 0.5* R - 0.4187 *G - 0.0813 *B#+128
return [Cb,Cr,Y]
def ycbcr2rgb3(dtuple):
Cb,Cr,Y=dtuple
R = Y + 1.402 *(Cr)#-128)
G = Y - 0.34414 *(Cb) - 0.71414 *(Cr)
B = Y + 1.772 *(Cb)
return np.clip([R,G,B],0.,255)
def rgb2YCbCr(img):
#print(img.shape)
img=np.array([[rgb2ycbcr3(y) for y in x] for x in img])
return img
def YCbCr2rgb(img):
img=np.array([[ycbcr2rgb3(y) for y in x] for x in img])
return img
def median_chrominance(img):
#HINT: ndimage.median_filter(image, size=10)
s=9
img=dstack((ndimage.median_filter(img[:,:,0],size=s),
ndimage.median_filter(img[:,:,1],size=s),
img[:,:,2]))
return img #!change accordingly
# This function applies a median chrominance filter on an RGB image
def filter_median_chrominance(img):
img = rgb2YCbCr(img)
img = median_chrominance(img)
img = YCbCr2rgb(img)
return img
lh_eb = demosaic_simple(img)
lh_f= filter_median_chrominance(lh_eb).astype(img.dtype)#/255.
lh_filtered =lh_f #white_balance_gray( lh_f)
show_images([lh_eb, lh_filtered], zoom=0.5, needs_encoding=True)
show_images([lh_eb[450:550, 50:150], lh_filtered[450:550, 50:150]], zoom=3.5, needs_encoding=True)
In [ ]: